/*=============================================================

PROGRAMMING TECHNIQUES FOR SUPERCOMPUTERS

Assignment Number 6

Task 3

=============================================================*/

#include <stdio.h>
#include <iostream>
#include "timing.hh"
#include <math.h>
#include <fstream>
#include <sys/time.h>
#include <malloc.h>
//#include <xmmintrin.h>
#include <nmmintrin.h>
#define Align 1

double timeStamp();
typedef double real;


int main()
{
  
  double wcTimeStart,wcTimeEnd,cpuTimeStart,cpuTimeEnd,wcTime, cpTime;	
	
  // Plot A for Vector Triad
  std::ofstream fileouta("SIMD.dat");
  fileouta << "# Plot for Task (a)" << std::endl;
  fileouta << "# N\t MFLOP/sec" << std::endl;
  
  std::ofstream fileoutb("noSIMD.dat");
  fileoutb << "# Plot for Task (b)" << std::endl;
  fileoutb << "# N\t MFLOP/sec" << std::endl;

  int adjustL=0;
  __m128d bi;
  __m128d ci;
  __m128d di;
  __m128d mult;
  __m128d rhs;
  for (int k =3; k<21; ++k)
  {
      int N = int(pow(2.5,k));
      //std::cout<<"Aligned\n";
      real *a = (real*) _mm_malloc(N*sizeof(real), 16);
      real *b = (real*) _mm_malloc(N*sizeof(real), 16);
      real *c = (real*) _mm_malloc(N*sizeof(real), 16);
      real *d = (real*) _mm_malloc(N*sizeof(real), 16);
      (N%2!=0)?adjustL=1:adjustL=0;
      /*for (int i = 0; i<N-adjustL; i+=2)
      {	
	 __m128d toAdd=_mm_set_pd(i+1.0,i+0.0);
	  _mm_store_pd(&b[i],toAdd);
	  _mm_store_pd(&c[i],toAdd);
	  _mm_store_pd(&d[i],toAdd);
      }
      if(adjustL!=0)
      {
	b[N-1]=N-1;
	c[N-1]=N-1;
	d[N-1]=N-1;
	}*/
      for (int i = 0; i<N; ++i)
      {	
	b[i] = i*2.0;
	c[i] = i*3.0;
	d[i] = i*1.0;
      }
      int repeat =1;
      real runtime =0.;
      for (; runtime <0.1; repeat*=2)
	{	
	  timing(wcTimeStart,cpuTimeStart);
	  for (int r=0; r<repeat; ++r)
	    {
	      //#pragma vector always
	      //#pragma simd
	      for (int i = 0; i<N-adjustL; i+=2)
		{
		   bi=_mm_load_pd(&b[i]);
		   ci=_mm_load_pd(&c[i]);
		   di=_mm_load_pd(&d[i]);
		   mult =_mm_mul_pd(ci,di);
		   rhs=_mm_add_pd(mult,bi);
		  _mm_store_pd(&a[i],rhs);
		}
	      if(adjustL!=0)
	      {
		a[N-1]=b[N-1]+c[N-1]*b[N-1];
	      }
	      if (a[N-1] <0.) std::cout << "Dummy Loop Encountered"<< std::endl;
	    }
	  timing(wcTimeEnd,cpuTimeEnd);	
	  runtime=wcTimeEnd-wcTimeStart;				
	}
      repeat/=2;
      long noFlops = N*2*repeat;
      real mFlops = noFlops/runtime/1000000;
      fileouta << N <<"\t" << mFlops << std::endl;
      _mm_free(a);
      _mm_free(b);
      _mm_free(c);
      _mm_free(d);
      
      double *a1 = new double[N];
      double *b1 = new double[N];
      double *c1 = new double[N];
      double *d1 = new double[N];
      for (int i = 0; i<N; ++i)
      {	
	b1[i] = i*2.0;
	c1[i] = i*3.0;
	d1[i] = i*1.0;
      }
      repeat =1;
      runtime =0.;
      for (; runtime <0.1; repeat*=2)
	{	
	  timing(wcTimeStart,cpuTimeStart);
	  for (int r=0; r<repeat; r++)
	    {
	      for (int i = 0; i<N; i++)
	      {
		a1[i] = b1[i] + c1[i]*d1[i];
	      }
	      if (a1[N-1] <0.) std::cout << "Dummy Loop Encountered"<< std::endl;
	    }
	  timing(wcTimeEnd,cpuTimeEnd);	
	  runtime=wcTimeEnd-wcTimeStart;				
	}
      repeat/=2;
      noFlops = N*2*repeat;
      mFlops = noFlops/runtime/1000000;
      fileoutb << N <<"\t" << mFlops << std::endl;
      delete [] a1;
      delete [] b1;
      delete [] c1;
      delete [] d1;
  }
  fileouta.close();
  fileoutb.close();
  return 0;
}

double timeStamp()
{
    struct timeval tp;
    gettimeofday(&tp, NULL);
    return((double)(tp.tv_sec + tp.tv_usec/1000000.0));
}
